4  Prognosemodel Retentie na 1 jaar

ITD | B Communication and Multimedia Design (CMD) - voltijd - versie 1.0

4.1 Methode, data en analyse

4.1.1 Toelichting op de methode

Voor de ontwikkeling van prognosemodellen gebruiken we de aanpak van Tidymodels. Tidymodels is een framework voor het bouwen van een prognosemodel. Hiermee verzekeren we ons van een systematische, herhaalbare en schaalbare aanpak. De code in deze pagina wordt getoond als deze hoort bij de ontwikkeling van een model. Overige code, voor bijvoorbeeld algemene tabellen of grafieken, kan zichtbaar gemaakt worden.

4.1.2 Toelichting op de data

De basis voor deze analyse is studiedata van De Haagse Hogeschool (De HHs), verrijkt door het lectoraat LTA. De data bevat informatie over de inschrijvingen van studenten in het eerste jaar van de opleiding:

  1. Demografische kenmerken: geslacht, leeftijd, reistijd en SES totaalscore.
  2. Vooropleidingskenmerken: toelaatgevende vooropleiding, studiekeuzeprofiel en gemiddeld eindcijfer in de vooropleiding.
  3. Aanmeldingskenmerken: aansluiting (direct na diploma, tussenjaar, switch), dag van aanmelding, aantal parallelle studies aan De HHs en collegejaar.

De variabelen die we hebben geselecteerd voor analyse is op basis van standaardgegevens van een popoluatie (leeftijd en geslacht), eerdere analsyes op voorspelkracht voor retentie na 1 jaar (Bakker, 2022), of omdat een kenmerk volgens Europese normen geldt als een sensitief kenmerk voor discriminatie (European Union Agency for Fundamental Rights, 2018).

Toon code
# Read the data dictionary
dfData_dictionary <- Get_Data_Dictionary()

# Show the data dictionary
Get_tblData_Dictionary(dfData_dictionary)
Tabel 4.1: Variabelen en mogelijke waarden

Variabele

Toelichting

Mogelijke waarden

Aanmelding

De dag van de aanmelding voor de studie gerekend vanaf 1 september

-366 - 366

Aansluiting

De manier waarop een student instroomt in de opleiding

2e studie, Direct, Na CD, Onbekend, Overig, Switch Extern, Switch Intern, Tussenjaar

APCG

Of de student ten tijde van het behalen van de toelaatgevende vooropleiding in een armoedeprobleemcumulatiegebied woonde

Ja, Nee, Onbekend

Cijfer_CE_Engels

Het gemiddelde cijfer Engels van het centaal examen van de middelbare school

0-10

Cijfer_CE_Natuurkunde

Het gemiddelde cijfer voor Natuurkunde van het centaal examen van de middelbare school

0-10

Cijfer_CE_Nederlands

Het gemiddelde cijfer Nederlands van het centaal examen van de middelbare school

0-10

Cijfer_CE_VO

Het gemiddelde cijfer voor het centaal examen van de middelbare school

0-10

Cijfer_CE_Wiskunde

Het gemiddelde cijfer voor Wiskunde van het centaal examen van de middelbare school

0-10

Cijfer_SE_VO

Het gemiddelde cijfer voor het schoolexamen van de middelbare school

0-10

Collegejaar

Het collegejaar van de inschrijving

2012-2022

Dubbele_studie

Gegeven of de student meer dan 1 studie volgt

Ja, Nee

Geslacht

Het geslacht van de student

M, V

ID

Uniek nummer per student

Leeftijd

De leeftijd van de student op 1 oktober

16-100

Ranking

De positie die de student heeft in de selectie voor de opleiding (alleen bij de B Huidtherapie)

0-350

Reistijd

De reistijd van de student naar de vestiging van de opleiding op een maandag om 9 uur met het openbaar vervoer vanaf de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding in minuten

0-120

Retentie

Of de student doorstudeert na het eerste studiejaar

Ja, Nee

SES_Arbeid

De sociaaleconomische status score op arbeid van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding

-1 - 1

SES_Totaal

De sociaaleconomische status score totaal van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding

-1 - 1

SES_Welvaart

De sociaaleconomische status score op welvaart van het CBS op basis van de postcode waarop de student woonde ten tijde van het behalen van de toelaatgevende vooropleiding

-1 - 1

Studiekeuzeprofiel

Het studiekeuzeprofiel in het voortgezet onderwijs (havo of vwo) of profiel in het mbo

AHO, ALG, CERT, CM, EA, EM, EM&CM, HB, HO, ICT, MedV, NG, NT, NT&NG, Onbekend, TP, TR, TSL, VNL, VS, ZW

Vooropleiding

De toelaatgevende vooropleiding

BD, CD, HAVO, HO, MBO, Overig, VWO

4.1.3 Toelichting op de analyse

We toetsen in deze analyse Retentie na 1 jaar, voortaan Retentie genoemd.

  • Retentie is gedefinieerd als ingeschreven staan in dezelfde opleiding in een aansluitend collegejaar. Een wisseling van opleidingsvorm binnen de opleiding, bijvoorbeeld van voltijd in jaar 1 naar duaal in jaar 2, geldt ook als retentie.
  • Uitval is het tegenovergestelde van retentie: niet ingeschreven staan in dezelfde opleiding in een aansluitend collegejaar.
  • Een wisseling van opleidingsvorm binnen de opleiding, bijvoorbeeld van voltijd in jaar 1 naar duaal in jaar 2, geldt niet als uitval.

4.2 Voorbereidingen

4.2.1 Laad de data

We laden een subset in van historische data specifiek voor:

Opleiding: ITD | B Communication and Multimedia Design (Synth) (CMD), voltijd, eerstejaars - Retentie na 1 jaar

Toon code
# Read the data for this study progamme
if(params$use_synthetic_data) {
  dfOpleiding_inschrijvingen_base <- Get_Studyprogram_Enrollments_Synthetic(
    studytrack = params$opleiding,
    studyform = toupper(params$opleidingsvorm_afkorting)) |> 
    mutate(
      ID = as.character(ID),
      INS_Opleiding = as.character(INS_Opleiding),
      INS_Opleidingsvorm = as.character(INS_Opleidingsvorm)
    ) |> 
    mutate(
      INS_Student_UUID_opleiding_vorm = paste(ID, INS_Opleiding, INS_Opleidingsvorm, sep = "_"),
      INS_Opleidingsnaam_huidig = paste(INS_Opleidingsnaam_huidig, "(Synth.)", sep = " ")
    )
} else {
  dfOpleiding_inschrijvingen_base <- get_lta_studyprogram_enrollments_pin(
    board = "HHs/Inschrijvingen",
    faculty = params$faculteit,
    studyprogram = params$opleidingsnaam,
    studytrack = params$opleiding,
    studyform = toupper(params$opleidingsvorm),
    range = "eerstejaars")
}

# Adjust dfOpleiding_inschrijvingen_base
dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |> 
  
  # Rearrange the levels
  mutate(across(all_of(names(lLevels_formal)), ~ factor(.x, 
                                                   levels = lLevels_formal[[cur_column()]]))) |> 

  # Create a simple success variable
  Mutate_Retentie(sSucces_model) |>
  
  # Convert the success variable into a factor
  mutate(SUC_Retentie = as.factor(SUC_Retentie)) |> 

  ## Special possibly based on the propaedeutic diploma
  # Filter_Propedeusediploma(sPropedeusediploma) |>

  # Make the Dual Study variable a Yes/No variable
  Mutate_Dubbele_studie() |>  

  # Remove INS_Aantal_inschrijvingen
  select(-INS_Aantal_inschrijvingen) 

  ## Adjust the levels of sensitive variables
  for(i in lSentitive_formal_variables){
    dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |>
      Mutate_Levels(
        i,
        list(lLevels_formal[[i]])
      )
  }
  
# B Huidtherapie: Filter on only students with a grade number (selection)
if(opleiding == "HDT") {
  dfOpleiding_inschrijvingen_base <- dfOpleiding_inschrijvingen_base |> 
    filter(!is.na(RNK_Rangnummer)) 
} 

4.2.2 Selecteer en inspecteer de data

We selecteren eerst de relevante variabelen. We verwijderen daarbij variabelen die maar één waarde hebben, omdat die geen voorspellende waarde kunnen hebben.

Toon code
lSelect <- Get_lSelect(dfVariables, "VAR_Formal_variable")

# B Huidtherapie: add the variable RNK_Rangnummer unless it is HDT
if(opleiding == "HDT") {
  lSelect <- c(lSelect, "RNK_Rangnummer")
}

# Create a mapping with formal and simple maes
lName_mapping <- setNames(dfVariables$VAR_Simple_variable, dfVariables$VAR_Formal_variable)

# Create a subset
dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen_base |>
  
  # Select the relevant variables
  select(all_of(lSelect)) |>
  
  # Rename variables for more readable names
  rename_with(~ lName_mapping[.x], .cols = everything()) |> 
  
  # Adjust CBS_APCG_tf to a factor
  Mutate_APCG() |>

  # Indicate where missing numbers are in VO
  Mutate_Cijfers_VO() |>
  
  # Remove variables, where there is only 1 value
  select(where(~ n_distinct(.) > 1)) |>
  
  # Sort
  arrange(Collegejaar, ID)

dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
 ltabase::sort_distinct()

## Remove the base dataset
#rm(dfOpleiding_inschrijvingen_base)

We inspecteren de variabelen in een samenvatting in relatie tot retentie en corrigeren daarbij voor multiple testing; de gecorrigeerde significantie-waarden staan vermeld als q-value.

Toon code
# Create a summary of the data
dfSummary <- dfOpleiding_inschrijvingen |>
  
  # Remove columns not relevant to the analysis
  select(-c(ID, Collegejaar)) |> 
  
  # Adjust the labels of Retention from True to Ja, and from False to Nee
  mutate(Retentie = fct_recode(Retentie, "Nee" = "FALSE", "Ja" = "TRUE")) |>
  
  # Adjust the order of the labels of Retentie
  mutate(Retentie = fct_relevel(Retentie, "Ja", "Nee")) |> 
  
  # Adjust Studiekeuzeprofiel: if NA, then unknown
  mutate(Studiekeuzeprofiel = coalesce(Studiekeuzeprofiel, "Onbekend")) |> 
  
  # Factor all character variables
  mutate(across(where(is.character), as.factor))

# Create a tbl_df of the summary
tbl_dfSummary <- Get_tblSummary(dfSummary) 

tbl_dfSummary
Tabel 4.2: Variabelen in relatie tot de uitkomstmaat: Retentie na 1 jaar

Retentie

Variabele

Ja, N=1002 (62%)1

Nee, N=611 (38%)1

p-value2

q-value3

Totaal, N = 16131

Aanmelding

135,80 (62,44)

112,36 (68,57)

<0,001

<0,001***

126,92 (65,80)

Aansluiting

<0,001

0,012***

Direct

455 (60%)

301 (40%)

756 (100%)

Tussenjaar

125 (71%)

51 (29%)

176 (100%)

Switch intern

148 (73%)

54 (27%)

202 (100%)

Switch extern

254 (57%)

191 (43%)

445 (100%)

2e Studie

6 (40%)

9 (60%)

15 (100%)

Na CD

14 (74%)

5 (26%)

19 (100%)

Overig

0 (NA%)

0 (NA%)

0 (NA%)

Onbekend

0 (NA%)

0 (NA%)

0 (NA%)

APCG

<0,001

0,013***

Ja

224 (63%)

133 (37%)

357 (100%)

Nee

721 (64%)

411 (36%)

1.132 (100%)

Onbekend

57 (46%)

67 (54%)

124 (100%)

Cijfer_CE_Engels

6,93 (1,18)

7,07 (1,10)

0,16

>0,99

6,98 (1,15)

Cijfer_CE_Engels_missing

0,70

>0,99

Ja

533 (62%)

331 (38%)

864 (100%)

Nee

469 (63%)

280 (37%)

749 (100%)

Cijfer_CE_Natuurkunde

6,18 (0,96)

6,31 (0,95)

0,21

>0,99

6,23 (0,95)

Cijfer_CE_Natuurkunde_missing

0,68

>0,99

Ja

867 (62%)

533 (38%)

1.400 (100%)

Nee

135 (63%)

78 (37%)

213 (100%)

Cijfer_CE_Nederlands

5,94 (0,85)

5,96 (0,91)

0,86

>0,99

5,95 (0,87)

Cijfer_CE_Nederlands_missing

0,61

>0,99

Ja

533 (62%)

333 (38%)

866 (100%)

Nee

469 (63%)

278 (37%)

747 (100%)

Cijfer_CE_VO

6,46 (0,37)

6,35 (0,34)

<0,001

<0,001***

6,42 (0,36)

Cijfer_CE_VO_missing

0,59

>0,99

Ja

488 (61%)

306 (39%)

794 (100%)

Nee

514 (63%)

305 (37%)

819 (100%)

Cijfer_CE_Wiskunde

6,35 (1,12)

6,28 (1,07)

0,65

>0,99

6,33 (1,10)

Cijfer_CE_Wiskunde_missing

0,64

>0,99

Ja

544 (62%)

339 (38%)

883 (100%)

Nee

458 (63%)

272 (37%)

730 (100%)

Cijfer_SE_VO

6,46 (0,38)

6,35 (0,36)

<0,001

<0,001***

6,42 (0,38)

Cijfer_SE_VO_missing

0,59

>0,99

Ja

501 (61%)

314 (39%)

815 (100%)

Nee

501 (63%)

297 (37%)

798 (100%)

Dubbele_studie

<0,001

<0,001***

Ja

12 (26%)

34 (74%)

46 (100%)

Nee

990 (63%)

577 (37%)

1.567 (100%)

Geslacht

<0,001

<0,001***

M

521 (57%)

394 (43%)

915 (100%)

V

481 (69%)

217 (31%)

698 (100%)

Leeftijd

20,01 (2,26)

20,19 (2,61)

0,48

>0,99

20,08 (2,40)

Reistijd

37,01 (18,38)

38,84 (22,79)

0,79

>0,99

37,70 (20,17)

SES_Arbeid

0,01 (0,08)

0,01 (0,08)

0,57

>0,99

0,01 (0,08)

SES_Totaal

0,02 (0,28)

0,01 (0,27)

0,45

>0,99

0,02 (0,28)

SES_Welvaart

0,01 (0,13)

0,01 (0,13)

0,36

>0,99

0,01 (0,13)

Studiekeuzeprofiel

0,044

>0,99*

AHO

2 (67%)

1 (33%)

3 (100%)

ALG

0 (0%)

1 (100%)

1 (100%)

CERT

1 (100%)

0 (0%)

1 (100%)

CM

78 (72%)

30 (28%)

108 (100%)

EA

36 (67%)

18 (33%)

54 (100%)

EM

118 (56%)

93 (44%)

211 (100%)

EM&CM

152 (70%)

66 (30%)

218 (100%)

HB

8 (62%)

5 (38%)

13 (100%)

HO

27 (57%)

20 (43%)

47 (100%)

ICT

82 (64%)

47 (36%)

129 (100%)

MedV

132 (65%)

72 (35%)

204 (100%)

NG

115 (63%)

68 (37%)

183 (100%)

NT

56 (60%)

38 (40%)

94 (100%)

NT&NG

62 (60%)

41 (40%)

103 (100%)

Onbekend

90 (52%)

83 (48%)

173 (100%)

TP

2 (40%)

3 (60%)

5 (100%)

TR

5 (56%)

4 (44%)

9 (100%)

TSL

4 (80%)

1 (20%)

5 (100%)

VNL

8 (80%)

2 (20%)

10 (100%)

VS

1 (33%)

2 (67%)

3 (100%)

ZW

23 (59%)

16 (41%)

39 (100%)

Vooropleiding

0,012

0,29*

MBO

330 (63%)

192 (37%)

522 (100%)

HAVO

548 (64%)

312 (36%)

860 (100%)

VWO

34 (59%)

24 (41%)

58 (100%)

BD

42 (46%)

50 (54%)

92 (100%)

CD

21 (70%)

9 (30%)

30 (100%)

HO

27 (53%)

24 (47%)

51 (100%)

Overig

0 (NA%)

0 (NA%)

0 (NA%)

Onbekend

0 (NA%)

0 (NA%)

0 (NA%)

1Mean (SD); n (%)

2Wilcoxon rank sum test; Fisher's Exact Test for Count Data with simulated p-value
(based on 2000 replicates); Pearson's Chi-squared test

3*p<0.05; **p<0.01; ***p<0.001

Daarnaast inspecteren we de kwaliteit van de data op missende waarden.

Toon code
# Load dlookr (temporary to avoid conflicts)
suppressMessages(library(dlookr))

# Show a summary of the data, sorted by missing values
diagnose(dfOpleiding_inschrijvingen) |> 
  mutate(missing_percent = round(missing_percent, 2),
         unique_rate = round(missing_percent, 2)) |>
  arrange(desc(missing_percent)) |>
  knitr::kable(col.names = c("Variabelen",
                           "Type",
                           "# Missende waarden",
                           "% Missende waarden",
                           "# Unieke waarden",
                           "Ratio unieke waarden"))
# Detach dlookr
detach("package:dlookr", unload = TRUE)
Tabel 4.3: Kwaliteit van de data voor bewerkingen (gesorteerd op missende waarden)
Variabelen Type # Missende waarden % Missende waarden # Unieke waarden Ratio unieke waarden
Cijfer_CE_Natuurkunde numeric 1400 86.79 42 86.79
Cijfer_CE_Wiskunde numeric 883 54.74 56 54.74
Cijfer_CE_Nederlands numeric 866 53.69 45 53.69
Cijfer_CE_Engels numeric 864 53.56 57 53.56
Cijfer_SE_VO numeric 815 50.53 23 50.53
Cijfer_CE_VO numeric 794 49.23 23 49.23
Studiekeuzeprofiel factor 173 10.73 21 10.73
SES_Welvaart numeric 126 7.81 355 7.81
SES_Arbeid numeric 125 7.75 261 7.75
SES_Totaal numeric 125 7.75 459 7.75
Reistijd numeric 20 1.24 362 1.24
Aanmelding numeric 0 0.00 266 0.00
Aansluiting factor 0 0.00 6 0.00
APCG character 0 0.00 3 0.00
Cijfer_CE_Engels_missing character 0 0.00 2 0.00
Cijfer_CE_Natuurkunde_missing character 0 0.00 2 0.00
Cijfer_CE_Nederlands_missing character 0 0.00 2 0.00
Cijfer_CE_VO_missing character 0 0.00 2 0.00
Cijfer_CE_Wiskunde_missing character 0 0.00 2 0.00
Cijfer_SE_VO_missing character 0 0.00 2 0.00
Collegejaar numeric 0 0.00 11 0.00
Dubbele_studie character 0 0.00 2 0.00
Geslacht factor 0 0.00 2 0.00
ID character 0 0.00 1613 0.00
Leeftijd numeric 0 0.00 17 0.00
Retentie factor 0 0.00 2 0.00
Vooropleiding factor 0 0.00 6 0.00

4.2.3 Bewerk de data

  • Uit de eerste diagnose blijkt dat niet alle variabelen goed genoeg zijn voor het bouwen van een prognosemodel: er zijn missende waarden en niet alle veldtypes zijn geschikt.
  • Om bias te voorkomen verwijderen we geen rijen met missende waarden, maar vullen die op (imputatie). We bewerken de data zo dat alle missende waarden worden opgevuld: bij numerieke waarden met het gemiddelde en bij categorische variabelen met ‘Onbekend’.
  • We passen het type van sommige variabelen aan, zodat ze in het model gebruikt kunnen worden: tekstvelden zetten we om naar factor (een categorische variabele); logische variabelen (Ja/Nee) zetten we om naar een numerieke variabele (1/0).
  • De uitkomstvariabele, Retentie, leiden we af van de variabele SUC_Uitval_aantal_jaar_LTA. Als de waarde daar 1 is, is de student na 1 jaar uitgevallen, 2 na 2 jaar, etc. Zolang de waarde daar 0 is, is de student niet uitgevallen.
  • Een fictief studentnummer (INS_Student_UUID_opleiding_vorm) gebruiken we, zodat we – als er afwijkende resultaten zijn – de dataset gericht kunnen onderzoeken als dat nodig is.
Toon code
# Edit the data
dfOpleiding_inschrijvingen <- dfOpleiding_inschrijvingen |> 
  
  # Imputate all numeric variables with the mean
  mutate(across(where(is.numeric), ~ ifelse(
    is.na(.x),
    mean(.x, na.rm = T),
    .x
  )) ) |>
  
  # Convert character variables to factor
  mutate(across(where(is.character), as.factor)) |> 
  
  # Convert logical variables to 0 or 1
  mutate(across(where(is.logical), as.integer)) |>
  
  # Fill in factors missing values with “Unknown”
  mutate(across(where(is.factor), ~ suppressWarnings(
    fct_explicit_na(.x, na_level = "Onbekend")
  ))) |> 
  
  # Rearrange the columns so that Retentie is in front
  select(Retentie, everything()) 

## View the data
# glimpse(dfOpleiding_inschrijvingen) 

# Load dlookr (temporary to avoid conflicts)
suppressMessages(library(dlookr))

# Diagnose the data
diagnose(dfOpleiding_inschrijvingen) |> 
  mutate(missing_percent = round(missing_percent, 2),
         unique_rate = round(unique_rate, 2)) |>
  knitr::kable(col.names = c("Variabelen",
                           "Type",
                           "# Missende waarden",
                           "% Missende waarden",
                           "# Unieke waarden",
                           "Ratio unieke waarden"))
# Detach dlookr
detach("package:dlookr", unload = TRUE)
Tabel 4.4: Kwaliteit van de data na bewerkingen (gesorteerd op missende waarden)
Variabelen Type # Missende waarden % Missende waarden # Unieke waarden Ratio unieke waarden
Retentie factor 0 0 2 0.00
Aanmelding numeric 0 0 266 0.16
Aansluiting factor 0 0 6 0.00
APCG factor 0 0 3 0.00
Cijfer_CE_Engels numeric 0 0 57 0.04
Cijfer_CE_Engels_missing factor 0 0 2 0.00
Cijfer_CE_Natuurkunde numeric 0 0 42 0.03
Cijfer_CE_Natuurkunde_missing factor 0 0 2 0.00
Cijfer_CE_Nederlands numeric 0 0 45 0.03
Cijfer_CE_Nederlands_missing factor 0 0 2 0.00
Cijfer_CE_VO numeric 0 0 23 0.01
Cijfer_CE_VO_missing factor 0 0 2 0.00
Cijfer_CE_Wiskunde numeric 0 0 56 0.03
Cijfer_CE_Wiskunde_missing factor 0 0 2 0.00
Cijfer_SE_VO numeric 0 0 23 0.01
Cijfer_SE_VO_missing factor 0 0 2 0.00
Collegejaar numeric 0 0 11 0.01
Dubbele_studie factor 0 0 2 0.00
Geslacht factor 0 0 2 0.00
ID factor 0 0 1613 1.00
Leeftijd numeric 0 0 17 0.01
Reistijd numeric 0 0 362 0.22
SES_Arbeid numeric 0 0 261 0.16
SES_Totaal numeric 0 0 459 0.28
SES_Welvaart numeric 0 0 355 0.22
Studiekeuzeprofiel factor 0 0 21 0.01
Vooropleiding factor 0 0 6 0.00

4.2.4 Inspecteer de onderlinge correlaties

Het is verstandig om voorafgaand aan het bouwen van een model te kijken naar de onderlinge correlaties tussen numerieke variabelen. Dit geeft inzicht in de data en kan helpen bij het maken van keuzes voor het model of de duiding van de uitkomsten.

We selecteren uitsluitend numerieke waarden en variabelen die een standaard deviatie hebben die groter is dan 0. We clusteren de data op basis van 4 clusters. De correlatie van de diagonaal verbergen we, aangezien deze altijd 1 is.

Clusters in correlaties kunnen per opleiding verschillen. Onderzoek in de correlaties hoe sterk deze zijn, welke clusters gevormd worden en of deze – vanuit de context van de opleiding, faculteit of onderwijsinstelling – logisch zijn. Verwerk de inzichten eventueel in een oplegger. Gangbare clusters zijn:

  • Cijfers: Een correlatie tussen Cijfer SO en Cijfer VO: dit geeft aan dat het schoolexamen hetzelfde meet als het centraal examen. Exacte vakken en talen kunnen met elkaar correleren.
  • Reistijd en SES: Deze zijn gecorreleerd aangezien SES scores zijn afgeleid van wijken. Wijken hebben een specifieke reistijd en positionering ten opzichte van de vestigingslocaties van De HHs kennen.
  • Leeftijd is vaak niet gecorreleerd met andere variabelen.
Toon code
# Create a plot of the intercorrelations in numerical variables
# Remove columns with a standard deviation of 0
dfCorrelation <- dfOpleiding_inschrijvingen |> 
  select(-Collegejaar) |>
  select(where(is.numeric)) |> 
  select_if(~ sd(.) != 0) |>
  cor()

dfCorrelation |>  
corrplot::corrplot(
  order = 'hclust', 
  addrect = 4,
  method = "number",
  hclust.method = "complete",
  tl.cex = 0.8,       
  tl.col = "black",
  diag = FALSE)
Figuur 4.1: Correlatiematrix
Toon code
# Apply hierarchical clustering
dist   <- as.dist(1 - dfCorrelation)      # Create a distance matrix
hclust <- hclust(dist, method = "complete")  # Hierarchical clustering

# Plot the dendrogram
# plot(hclust, cex = 0.8)

# Create a clustering
dfClusters <- cutree(hclust, k = 4)

# Show clustering
# dfClusters

4.2.5 Bouw de trainingset, validatieset en testset

  • De data is nu geschikt om een prognosemodel mee te bouwen.
  • Om het model te bouwen, testen en valideren, splitsen we de data in drie delen van 60%, 20% en 20%. We doen dit op zo’n manier, dat elk deel ongeveer een gelijk aantal studenten bevat dat doorstudeert (dus niet uitvalt).
  • We trainen het model op basis van 60% en valideren de modellen tijdens het trainen op de overige 20% (de validatieset).
  • De verdeling van de training- en validatieset muteren we 10x (10 folds) om te voorkomen dat het model te veel leert van de trainingset en daardoor slecht presteert op de validatieset.
  • Als het model klaar is, testen we het op de 20% studenten uit de testset. De testset blijft dus de gehele tijd ongemoeid, zodat we overfitting - een te goed model op bekende data, maar slechte presetaties (performance) op onbekende data - voorkomen.
  • Een willekeurig, maar vaststaand seed-getal voorkomt dat we bij elke run van het model c.q. deze code een net iets andere uitkomst krijgen.

We ontwikkelen in de verdere analyse eerst een aantal modellen en testen de prestaties daarvan op de trainingset en validiatieset. Vervolgens bepalen we welk model het beste presteert en passen dat model dan toe op de testset (de uiteindelijke fit). Vandaar dat het toetsen van de prestaties van de modellen meerdere keren behandeld wordt.

Toon code
knitr::include_graphics(here::here("R/images", "prognosemodel-dataset-lta-hhs.png"))
Figuur 4.2: Splitsing van de dataset in trainingset, validatieset en testset
Toon code
set.seed(0821)

# Split the data into 3 parts: 60%, 20% and 20%
splits      <- initial_validation_split(dfOpleiding_inschrijvingen,
                                        strata = Retentie,
                                        prop = c(0.6, 0.2))

# Create three sets: a training set, a test set and a validation set
dfRetentie_train      <- training(splits)
dfRetentie_test       <- testing(splits)
dfRetentie_validation <- validation_set(splits)

# Create a resample set based on 10 folds (default)
dfRetentie_resamples  <- vfold_cv(dfRetentie_train, strata = Retentie)
Tabel 4.5: Verhouding van de uitkomstvariabele in de training- en testset
Naam Retentie Aantal Proportie
Trainingset FALSE 366 37.8%
Trainingset TRUE 601 62.2%
Testset FALSE 123 38.0%
Testset TRUE 201 62.0%

4.3 Model I: Logistische Regressie

  • Het eerste model is een logistische regressie met penalized likelihood; we gebruiken de glmnet engine voor het bouwen van het model. Penalized likelihood is een techniek die helpt bij het voorkomen van overfitting: het voegt voor elke extra variabele een strafterm toe om eenvoudige modellen te belonen. Glmnet is een veelgebruikt package voor het bouwen van logistische regressiemodellen.
  • We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric. De ROC-curve (Receiver Operating Characteristic) is een grafiek die de prestaties van een classificatiemodel afbeeldt door de verhouding tussen de true positives (sensitiviteit) en de false positives (aspecficiteit = 1-specificiteit) te plotten bij verschillende drempelwaarden. De oppervlakte onder deze curve, bekend als de AUC (Area Under the Curve), kwantificeert het onderscheidingsvermogen van het model; een AUC van 1 duidt op een perfect onderscheidend vermogen, terwijl een AUC van 0,5 wijst op een model zonder onderscheidend vermogen.

4.3.1 Maak het model

Eerst bouwen we het model.

# Build the model: logistic regression
lr_mod <- 
  logistic_reg(penalty = tune(), mixture = 1) |> 
  set_engine("glmnet")

4.3.2 Maak de recipe

Vervolgens zetten we meerdere stappen in een ‘recipe’:

  • We definiëren de student-ID als ID variabele. Daarmee krijgt deze variabele de rol van uniek rij-kenmerk.
  • We verwijderen vervolgens de oorspronkelijke student-ID en het collegejaar uit de data, omdat deze verder niet gebruikt moeten worden in het model.
  • We converteren factoren naar dummy variabelen: voor elke categorie wordt er een nieuwe logische variabele (Ja/Nee) aangemaakt.
  • We verwijderen variabelen die geen waarde toevoegen: variabelen met uitsluitend nullen.
  • We normaliseren numerieke variabelen om ze met elkaar te kunnen vergelijken door ze te centreren en schalen: het transformeert numerieke gegevens zodat ze een standaard deviatie van één en een gemiddelde van nul hebben.
  • Sterk gecorreleerde waarden verwijderen we nu niet, omdat we later in de analyse de eventuele samenhang met andere variabelen in een prognosemodel nog willen kunnen visualiseren.
# Build the recipe: logistic regression
lr_recipe <- 
  recipe(Retentie ~ ., data = dfRetentie_train) |>  
  update_role(ID, new_role = "ID") |>           # Set the student ID as an ID variable
  step_rm(ID, Collegejaar) |>                   # Remove ID and college year from the model
  step_unknown(Studiekeuzeprofiel, 
               new_level = "Onbekend skp") |>   # Add unknown skp
  step_dummy(all_nominal_predictors()) |>       # Create dummy variables from categorical variables
  step_zv(all_predictors()) |>                  # Remove zero values
  step_normalize(all_numeric_predictors())      # Center and scale numeric variables
Toon code
# Show the recipe
tidy(lr_recipe) |> 
  knitr::kable(col.names = c("Nummer", 
                             "Operatie", 
                             "Type",
                             "Getraind",
                             "Sla over",
                             "ID"))
Tabel 4.6: Recipesteps voor logistische regressie
Nummer Operatie Type Getraind Sla over ID
1 step rm FALSE FALSE rm_t7fOp
2 step unknown FALSE FALSE unknown_KV8Lp
3 step dummy FALSE FALSE dummy_bUMsj
4 step zv FALSE FALSE zv_beaLb
5 step normalize FALSE FALSE normalize_3WTmw

De variabelen die nu nog overblijven zijn:

Tabel 4.7: Resterende variabelen voor logistische regressie na bewerkingen
Aanmelding APCG_Nee Studiekeuzeprofiel_HB
Cijfer_CE_Engels APCG_Onbekend Studiekeuzeprofiel_HO
Cijfer_CE_Natuurkunde Cijfer_CE_Engels_missing_Nee Studiekeuzeprofiel_ICT
Cijfer_CE_Nederlands Cijfer_CE_Natuurkunde_missing_Nee Studiekeuzeprofiel_MedV
Cijfer_CE_VO Cijfer_CE_Nederlands_missing_Nee Studiekeuzeprofiel_TP
Cijfer_CE_Wiskunde Cijfer_CE_VO_missing_Nee Studiekeuzeprofiel_TR
Cijfer_SE_VO Cijfer_CE_Wiskunde_missing_Nee Studiekeuzeprofiel_TSL
Leeftijd Cijfer_SE_VO_missing_Nee Studiekeuzeprofiel_VNL
Reistijd Dubbele_studie_Nee Studiekeuzeprofiel_VS
SES_Arbeid Geslacht_V Studiekeuzeprofiel_ZW
SES_Totaal Studiekeuzeprofiel_CM Studiekeuzeprofiel_Onbekend
SES_Welvaart Studiekeuzeprofiel_EM.CM Vooropleiding_HAVO
Retentie Studiekeuzeprofiel_NT Vooropleiding_VWO
Aansluiting_Tussenjaar Studiekeuzeprofiel_NG Vooropleiding_BD
Aansluiting_Switch.intern Studiekeuzeprofiel_NT.NG Vooropleiding_CD
Aansluiting_Switch.extern Studiekeuzeprofiel_AHO Vooropleiding_HO
Aansluiting_X2e.Studie Studiekeuzeprofiel_ALG
Aansluiting_Na.CD Studiekeuzeprofiel_EA

4.3.3 Maak de workflow

Voor de uitvoering bouwen we een workflow. Daaraan voegen we het model en de bewerkingen in de recipe toe.

# Create the workflow: logistic regression
lr_workflow <- 
  workflow() |>         # Create a workflow
  add_model(lr_mod) |>  # Add the model
  add_recipe(lr_recipe) # Add the recipe

# Show workflow
lr_workflow
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: logistic_reg()

── Preprocessor ────────────────────────────────────────────────────────────────
5 Recipe Steps

• step_rm()
• step_unknown()
• step_dummy()
• step_zv()
• step_normalize()

── Model ───────────────────────────────────────────────────────────────────────
Logistic Regression Model Specification (classification)

Main Arguments:
  penalty = tune()
  mixture = 1

Computational engine: glmnet 

4.3.4 Tune en train het model

Het model moet getuned worden. Dit houdt in dat we de beste parameters voor het model moeten vinden. We maken een grid met verschillende penalty waarden. Daarmee kunnen we vervolgens het beste model selecteren met de hoogste ROC/AUC. We plotten de resultaten van de tuning, zodat we hieruit het beste model kunnen kiezen.

# Create a grid: logistic regression
lr_reg_grid <- tibble(penalty = 10 ^ seq(-4, -1, length.out = 30))

# Train and tune the model: logistic regression
lr_res <- 
  lr_workflow |> 
  tune_grid(dfRetentie_validation,
            grid = lr_reg_grid,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
Toon code
# Plot the results + a red vertical line for the max AUC
lr_plot <- 
  lr_res |> 
  collect_metrics() |> 
  ggplot(aes(x = penalty, y = mean)) + 
  geom_point() + 
  geom_line() + 
  
  # Make the scale of the x-axis logarithmic
  scale_x_log10(labels = scales::label_number()) +
    theme(
      axis.title.x = element_text(margin = margin(t = 20))
    ) +
  
  # Define the title, subtitle and caption
  labs(
    caption = sCaption,
    x = "Area under the ROC Curve",
    y = "Penalty"
  )
  
  # Add LTA elements
  lr_plot <- Add_LTA_Theme_Elements(lr_plot, title_subtitle = FALSE)
  
# Find the penalty value with the max AUC
max_auc_penalty <- lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |> 
  pull(penalty)

# Add the red vertical line to lr_plot
lr_plot_plus <- lr_plot + 
  geom_vline(xintercept = max_auc_penalty, color = "red")

# Find a mean for the max AUC that is higher
max_auc_mean <- lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |> 
  pull(penalty)

# Print the final plot
lr_plot_plus
Figuur 4.3: Tuning resultaten logistische regressie

4.3.5 Kies het beste model

De prestaties van een model gevisualiseerd met behulp van een ROC curve. De sensitiviteit (True Positive Rate) en specificiteit (True Negative Rate) worden hierin tegenover elkaar uitgezet. De Area under the ROC Curve (AUC/ROC) geeft de prestaties van het model weer. Het model scoort beter naarmate de AUC/ROC dichter bij de 1 ligt, de linker bovenhoek. De linker bovenhoek houdt in dat alle prognoses exact overeenstemmen met de werkelijkheid. Een AUC/ROC van 0,5 betekent dat het model niet beter presteert dan een willekeurige voorspelling.

We gebruiken modellen met een zo hoog mogelijke Area under the ROC Curve (AUC/ROC) en een zo laag mogelijke penalty. Zo kunnen we uit de resultaten het beste model kiezen en visualiseren.

# Show the best model
top_models <-
  lr_res |> 
  show_best(metric = "roc_auc", n = 10) |> 
  mutate(mean = round(mean, 6)) |>
  arrange(penalty) 
Toon code
top_models|> 
  knitr::kable(col.names = c("Penalty", 
                             "Metriek", 
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))
Tabel 4.8: Model performance voor logistische regressie
Penalty Metriek Estimator Gemiddelde Aantal SE Configuratie
0.0017433 roc_auc binary 0.669959 1 NA Preprocessor1_Model13
0.0022122 roc_auc binary 0.671680 1 NA Preprocessor1_Model14
0.0028072 roc_auc binary 0.672992 1 NA Preprocessor1_Model15
0.0035622 roc_auc binary 0.674549 1 NA Preprocessor1_Model16
0.0045204 roc_auc binary 0.675041 1 NA Preprocessor1_Model17
0.0057362 roc_auc binary 0.676516 1 NA Preprocessor1_Model18
0.0072790 roc_auc binary 0.678566 1 NA Preprocessor1_Model19
0.0092367 roc_auc binary 0.678975 1 NA Preprocessor1_Model20
0.0117210 roc_auc binary 0.675492 1 NA Preprocessor1_Model21
0.0148735 roc_auc binary 0.668443 1 NA Preprocessor1_Model22
# Select the best model: logistic regression
lr_best <- 
  lr_res |> 
  collect_metrics() |> 
  filter(mean == max(mean)) |>
  slice(1) 
Toon code
lr_best|> 
  mutate(mean = round(mean, 6)) |>
  knitr::kable(col.names = c("Penalty", 
                             "Metriek", 
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))
Tabel 4.9: Hoogste model performance voor logistische regressie
Penalty Metriek Estimator Gemiddelde Aantal SE Configuratie
0.0092367 roc_auc binary 0.678975 1 NA Preprocessor1_Model20
# Collect the predictions and evaluate the model (AUC/ROC): logistic regression
lr_auc <- 
  lr_res |> 
  collect_predictions(parameters = lr_best) |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Logistisch Regressie")
Toon code
# Plot the ROC curve
Get_ROC_Plot(lr_auc, position = 1)
Figuur 4.4: ROC curve voor logistische regressie
Toon code
# Determine the AUC of the best model
lr_auc_highest   <-
  lr_res |>
  collect_predictions(parameters = lr_best) |> 
  roc_auc(Retentie, .pred_FALSE)

# Add model name and AUC dfModel_results
dfModel_results <- 
  dfModel_results |>
  add_row(model = "Logistic Regression", auc = lr_auc_highest$.estimate)

4.4 Model II: Tree-based ensemble

  • Het tweede model is een random forest: een ensemble van beslisbomen (decision trees). Het is een krachtig model dat goed om kan gaan met complexe data en veel variabelen.
  • We gebruiken de ranger engine voor het bouwen van het model.

4.4.1 Bepaal het aantal PC-cores

Omdat een random forest model veel berekeningen vereist, willen we daarvoor alle computerkracht gebruiken die beschikbaar is. Het aantal CPU’s (cores), wat verschilt per computer, bepaalt hoe snel het model getraind kan worden. We bepalen het aantal cores en gebruiken dat bij het bouwen van het model.

Toon code
# Determine the number of cores
cores <- parallel::detectCores()

4.4.2 Maak het model

We bouwen eerst het model. We gebruiken de rand_forest functie om het model te bouwen. We tunen de mtry en min_n parameters. De mtry parameter bepaalt het aantal variabelen dat per boom wordt gebruikt. De min_n parameter bepaalt het minimum aantal observaties dat in een blad van de boom moet zitten. De functie tune() is hier nog een placeholder om de beste waarden voor deze parameters - die we later bepalen - in te kunnen stellen. We gebruiken 1.000 bomen c.q. versies van het model.

# Build the model: random forest

rf_mod <- 
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) |> 
  set_engine("ranger", num.threads = cores) |> 
  set_mode("classification")

4.4.3 Maak de recipe

We maken een recipe voor het random forest model. We verwijderen de student ID en het collegejaar uit de data, omdat deze niet moet worden gebruikt in het model. Overige stappen zijn bij een random forest minder relevant in tegenstelling tot een regressiemodel.

# Create the recipe: random forest
rf_recipe <- 
  recipe(Retentie ~ ., data = dfRetentie_train) |> 
  step_unknown(Studiekeuzeprofiel, 
               new_level = "Onbekend skp") |>   # Add unknown skp
  step_rm(ID, Collegejaar)                      # Remove ID and Collegejaar from the model
Toon code
# Show the recipe
tidy(rf_recipe) |> 
  knitr::kable(col.names = c("Nummer", 
                             "Operatie", 
                             "Type",
                             "Getraind",
                             "Sla over",
                             "ID"))
Tabel 4.10: Recipesteps voor random forest
Nummer Operatie Type Getraind Sla over ID
1 step unknown FALSE FALSE unknown_QY996
2 step rm FALSE FALSE rm_Q9GeT

4.4.4 Maak de workflow

We voegen het model en de recipe toe aan de workflow voor dit model.

# Create the workflow: random forest
rf_workflow <- 
  workflow() |> 
  add_model(rf_mod) |> 
  add_recipe(rf_recipe)

# Show workflow
rf_workflow
══ Workflow ════════════════════════════════════════════════════════════════════
Preprocessor: Recipe
Model: rand_forest()

── Preprocessor ────────────────────────────────────────────────────────────────
2 Recipe Steps

• step_unknown()
• step_rm()

── Model ───────────────────────────────────────────────────────────────────────
Random Forest Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = 1000
  min_n = tune()

Engine-Specific Arguments:
  num.threads = cores

Computational engine: ranger 

4.4.5 Tune en train het model

We trainen en tunen het model in de workflow. We maken een grid met verschillende waarden voor de parameters mtry en min_n. We gebruiken de Area under the ROC Curve (AUC/ROC) als performance metric. Met de resultaten van de tuning kiezen we het beste model.

# Show the parameters that can be tuned
rf_mod
Random Forest Model Specification (classification)

Main Arguments:
  mtry = tune()
  trees = 1000
  min_n = tune()

Engine-Specific Arguments:
  num.threads = cores

Computational engine: ranger 
# Extract the parameters being tuned
extract_parameter_set_dials(rf_mod)
Collection of 2 parameters for tuning

 identifier  type    object
       mtry  mtry nparam[?]
      min_n min_n nparam[+]

Model parameters needing finalization:
   # Randomly Selected Predictors ('mtry')

See `?dials::finalize` or `?dials::update.parameters` for more information.
# Determine the seed
set.seed(2904)

# Build the grid: random forest
rf_res <- 
  rf_workflow |> 
  tune_grid(dfRetentie_validation,
            grid = 25,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
i Creating pre-processing data to finalize unknown parameter: mtry

4.4.6 Kies het beste model

We evalueren de beste modellen en maken een ROC curve om de performance van het model te visualiseren. Vervolgens vergelijken we de prestaties van de modellen en kiezen daaruit het beste model.

# Show the best models
rf_res |> 
  show_best(metric = "roc_auc", n = 15) |> 
  mutate(mean = round(mean, 6)) |>
  knitr::kable(col.names = c("Mtry", 
                             "Min. aantal", 
                             "Metriek",
                             "Estimator",
                             "Gemiddelde",
                             "Aantal",
                             "SE",
                             "Configuratie"))
Tabel 4.11: Model performance voor random forest
Mtry Min. aantal Metriek Estimator Gemiddelde Aantal SE Configuratie
4 26 roc_auc binary 0.698525 1 NA Preprocessor1_Model23
2 32 roc_auc binary 0.697869 1 NA Preprocessor1_Model03
7 35 roc_auc binary 0.697869 1 NA Preprocessor1_Model14
8 31 roc_auc binary 0.696393 1 NA Preprocessor1_Model24
22 38 roc_auc binary 0.695656 1 NA Preprocessor1_Model01
18 40 roc_auc binary 0.695533 1 NA Preprocessor1_Model16
5 36 roc_auc binary 0.695369 1 NA Preprocessor1_Model09
6 23 roc_auc binary 0.695164 1 NA Preprocessor1_Model17
3 8 roc_auc binary 0.693402 1 NA Preprocessor1_Model08
13 25 roc_auc binary 0.692172 1 NA Preprocessor1_Model25
16 34 roc_auc binary 0.691885 1 NA Preprocessor1_Model07
11 29 roc_auc binary 0.691721 1 NA Preprocessor1_Model02
9 12 roc_auc binary 0.691721 1 NA Preprocessor1_Model05
12 13 roc_auc binary 0.690697 1 NA Preprocessor1_Model19
11 21 roc_auc binary 0.690246 1 NA Preprocessor1_Model15
Toon code
# Plot the results
autoplot <- autoplot(rf_res) +
  theme_minimal() +
  labs(
    y = "roc/auc",
    caption = sCaption
  )
  
  # Add LTA elements
  autoplot <- Add_LTA_Theme_Elements(autoplot, title_subtitle = FALSE)
  
  print(autoplot)
Figuur 4.5: Model performance random forest
# Select the best model
rf_best <- 
  rf_res |> 
  select_best(metric = "roc_auc")
Toon code
rf_best|> 
  knitr::kable(col.names = c("Mtry", 
                             "Min. aantal", 
                             "Configuratie"))
Tabel 4.12: Hoogste model performance voor random forest
Mtry Min. aantal Configuratie
4 26 Preprocessor1_Model23
Toon code
# Collect the predictions
rf_res |> 
  collect_predictions() |> 
  head(10) |>
  mutate(.pred_FALSE = scales::percent(.pred_FALSE, accuracy = 0.1),
         .pred_TRUE = scales::percent(.pred_TRUE, accuracy = 0.1)) |>
  knitr::kable(col.names = c("% Voorsp. FALSE", 
                             "% Voorsp. TRUE", 
                             "ID",
                             "Rij",
                             "Mtry", 
                             "Min. aantal", 
                             "Retentie",
                             "Configuratie"))
Tabel 4.13: Predicties voor random forest
% Voorsp. FALSE % Voorsp. TRUE ID Rij Mtry Min. aantal Retentie Configuratie
52.6% 47.4% validation 968 22 38 TRUE Preprocessor1_Model01
42.9% 57.1% validation 969 22 38 TRUE Preprocessor1_Model01
40.4% 59.6% validation 970 22 38 TRUE Preprocessor1_Model01
39.5% 60.5% validation 971 22 38 TRUE Preprocessor1_Model01
30.4% 69.6% validation 972 22 38 TRUE Preprocessor1_Model01
49.0% 51.0% validation 973 22 38 TRUE Preprocessor1_Model01
21.2% 78.8% validation 974 22 38 FALSE Preprocessor1_Model01
48.2% 51.8% validation 975 22 38 TRUE Preprocessor1_Model01
50.1% 49.9% validation 976 22 38 FALSE Preprocessor1_Model01
27.4% 72.6% validation 977 22 38 TRUE Preprocessor1_Model01
Toon code
# Determine the AUC/ROC curve
rf_auc <- 
  rf_res |> 
  collect_predictions(parameters = rf_best) |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Random Forest")

# Plot the ROC curve
Get_ROC_Plot(rf_auc, position = 2)

# Determine the AUC of the best model
rf_auc_highest   <-
  rf_res |>
  collect_predictions(parameters = rf_best) |> 
  roc_auc(Retentie, .pred_FALSE)

# Add model name and AUC to dfModel_results
dfModel_results <- 
  dfModel_results |>
  add_row(model = "Random Forest", 
          auc = rf_auc_highest$.estimate)
Figuur 4.6: ROC curve voor random forest

4.5 De uiteindelijke fit

  • In de laatste stap van deze analyse maken we het model definitief.
  • We testen het model op de testset en evalueren het model met metrieken en de Variable Importance (VI). De VI kwantificeert de bijdrage van elke variabele aan de voorspellende kracht van een model. Het identificeert welke variabelen significant zijn voor de modelprestaties, wat essentieel is voor het interpreteren en optimaliseren van een model (Van der Laan, 2006). Methoden zoals de Shapley-waarde en permutation importance worden vaak toegepast om dit belang te meten. Op deze methoden komen we terug in het volgende hoofdstuk.

4.5.1 Combineer de AUC/ROC curves en kies het beste model

Eerst combineren we de AUC/ROC curves van de modellen om ze te vergelijken. We kiezen het beste model op basis van de hoogste AUC/ROC.

Toon code
# Combine the AUC/ROC curves to compare the models
Get_ROC_Plot(list(lr_auc, rf_auc))
Figuur 4.7: Gecombineerde ROC curves
Toon code
# Determine which of the models is best based on highest AUC/ROC
dfModel_results <- dfModel_results |>
  mutate(number = row_number()) |> 
  mutate(best = ifelse(auc == max(auc), TRUE, FALSE)) |> 
  arrange(number)

# Determine the best model
sBest_model     <- dfModel_results$model[dfModel_results$best == TRUE]
sBest_model_auc <- round(dfModel_results$auc[dfModel_results$best == TRUE], 4)

Het beste model is het Random Forest model met een AUC/ROC van 0.6985. Het Logistic Regression model heeft een AUC van 0.679. We ronden de analyse verder af met het Random Forest model op de validatieset.

4.5.2 Maak het finale model

We maken het finale model op basis van de beste parameters die we hebben gevonden. Door in de engine bij importance de impurity op te geven, wordt het beste random forest model gekozen om de data definitief mee te classificeren.

# Test the developed model on the test set
# Determine the optimal parameters

# Build the final models
last_lr_mod <-
    logistic_reg(penalty = lr_best$penalty,
                 mixture = 1) |>
    set_engine("glmnet") |>
    set_mode("classification")

last_rf_mod <-
    rand_forest(mtry = rf_best$mtry,
                min_n = rf_best$min_n,
                trees = 1000) |>
    set_engine("ranger", num.threads = cores, importance = "impurity") |>
    set_mode("classification")

4.5.3 Maak de workflow

We voegen het model toe aan de workflow en updaten de workflow met het finale model.

# Update the workflows
 last_lr_workflow <- 
    lr_workflow |> 
    update_model(last_lr_mod)

 last_rf_workflow <- 
    rf_workflow |> 
    update_model(last_rf_mod)

4.5.4 Fit het finale model

We voeren de finale fit uit. De functie last_fit past het model toe op de validatieset.

# Perform the final fit
set.seed(2904)

# Make a final fit for both models so we can save it for later use
last_fit_lr <- 
    last_lr_workflow |> 
    last_fit(splits)

last_fit_rf <- 
    last_rf_workflow |> 
    last_fit(splits)

lLast_fits <- list(last_fit_lr, last_fit_rf) |> 
  set_names(c("Logistic Regression", "Random Forest"))

# Determine which model is best
if(sBest_model == "Logistic Regression") {
  last_fit <- last_fit_lr
} else if(sBest_model == "Random Forest") {
  last_fit <- last_fit_rf
}

# Keep results, model results and associated data
sFittedmodels_outputpath <- Get_Model_Outputpath(mode = "last-fits")
saveRDS(lLast_fits, file = sFittedmodels_outputpath)

sModelresults_outputpath <- Get_Model_Outputpath(mode = "modelresults")
saveRDS(dfModel_results, file = sModelresults_outputpath)

sData_outputpath <- Get_Model_Outputpath(mode = "data")
saveRDS(dfOpleiding_inschrijvingen, file = sData_outputpath)

4.5.5 Evalueer het finale model: metrieken en variable importance

We evalueren het finale model nu grondiger op basis van 4 metrieken: 1) accuraatheid, 2) ROC/AUC en 3) de Brier score (de Mean Squared Error). Het is zinvol om accuraatheid, ROC/AUC en de Brier-score pas bij het finale model toe te passen, omdat dit efficiënter is en overfitting voorkomt. Zo combineren we een snelle modelselectie met een grondige evaluatie van het uiteindelijke model.

# Collect the metrics
last_fit |> 
  collect_metrics() |> 
  mutate(.estimate = round(.estimate, 4)) |>
  knitr::kable(col.names = c("Metriek", 
                             "Estimator",
                             "Estimate",
                             "Configuratie"))
Metriek Estimator Estimate Configuratie
accuracy binary 0.6420 Preprocessor1_Model1
roc_auc binary 0.6642 Preprocessor1_Model1
brier_class binary 0.2189 Preprocessor1_Model1

4.5.6 Plot de ROC curve

Tot slot visualiseren we de prestaties weer met een ROC curve van het beste model.

# Show the roc curve
auc_lf <- last_fit |> 
  collect_predictions() |> 
  roc_curve(Retentie, .pred_FALSE) |> 
  mutate(model = "Last fit")

Get_ROC_Plot(auc_lf, position = 3)
Figuur 4.8: ROC curve finale model

4.6 Conclusies

4.6.1 Het beste prognosemodel voor deze opleiding

Het beste prognosemodel blijkt het Random Forest model te zijn.

  • Van de prognosemodellen die we hebben ontwikkeld om retentie na 1 jaar te voorspellen, had het Random Forest model de hoogste AUC/ROC waarde (0.6985).

4.6.2 Mate van accuraatheid en lift

Een prognosemodel moet minimaal beter presteren dan een advanced-report om waarde op basis van accuraatheid toe te voegen. Het advanced-report neemt als basis de grootste klasse van de gemiddelde retentie na 1 jaar van de afgelopen jaren. Stel we zouden tegen alle studenten zeggen dat ze hun studie gaan halen, dan is de mate van accuratesse gelijk aan dit advanced-report. Dit advanced-report is dus altijd hoger dan de 50% lijn van de AUC/ROC curve, tenzij het advanced-report toevallig precies 50% is.

Toon code
knitr::include_graphics(here::here("R/images", "basemodel-lift.png"))
Figuur 4.9: Lift afhankelijk van advanced-report en accuraatheid

De mate van accuraatheid van het prognosemodel is vrij laag (64.2%).

  • advanced-report: 62.12% – Voor deze opleiding berekenen we het advanced-report als volgt. Van alle studenten studeerde 62.12% door; 100% - 62.12% = 37.88% studeerde niet door. De grootste klasse van deze twee, 62.12%, is daarmee de accuratesse van het advanced-report.
  • Accuratesse prognose: 64.2% – Het model voorspelt Retentie na 1 jaar met een accuratesse van 64.2%.
  • Lift: 2.08% – Het model scoort in de huidige opbouw met een verschil van 2.08% (de lift) iets beter dan de accuraatheid van het advanced-report.

4.6.3 Confusion Matrix

Toon code
# Determine the confusion matrix
confusion_matrix <- last_fit |>
  collect_predictions() |>
  conf_mat(truth = Retentie, estimate = .pred_class) 

dfConf_matrix <- as_tibble(confusion_matrix$table) |>
  rename(Werkelijkheid = Truth) |>
  mutate(Werkelijkheid = ifelse(Werkelijkheid == "TRUE", "Retentie", "Geen retentie"),
         Prediction    = ifelse(Prediction == "TRUE", "Retentie", "Geen retentie"))

pTP  <- Change_Number_Marks((dfConf_matrix$n[4]/sum(dfConf_matrix$n)*100),1)
pFP  <- Change_Number_Marks((dfConf_matrix$n[2]/sum(dfConf_matrix$n)*100),1)
pTN  <- Change_Number_Marks((dfConf_matrix$n[1]/sum(dfConf_matrix$n)*100),1)
pFN  <- Change_Number_Marks((dfConf_matrix$n[3]/sum(dfConf_matrix$n)*100),1)
pACC <- Change_Number_Marks(Last_fit_Accuracy,1)

De prestaties van het model kunnen we verder uitdrukken in een confusion matrix. Hierin zien we de voorspellingen van het model en de werkelijke uitkomsten. De matrix geeft inzicht in de mate van correcte en incorrecte voorspellingen. Ter illustratie werken we de matrix uit voor een voorspelling waarop een bindend studieadvies (BSA) gebaseerd zou kunnen zijn.

Toon code
knitr::include_graphics(here::here("R/images", "confusion-matrix-retention-lta-hhs.png"))
Figuur 4.10: Confusion matrix in relatie tot BSA

We passen de confusion matrix nu toe op het model dat als beste naar voren kwam. De accuraatheid van dit model is 64,2%. De accuraatheid van het model berekenen we door de som van de diagonaal te berekenen: het aandeel goed voorspelde uitkomsten, Retentie = Retentie (True Positive) en Geen retentie = Geen retentie (True Negative), af te zetten tegen het totaal aantal voorspellingen: 56,5% + 7,7% = 64,2%. (NB. De weergave in deze confusion matrix is diagonaal gespiegeld vergeleken met het voorbeeld.)

Toon code
confusion_plot <- plot_confusion_matrix(
    dfConf_matrix,
    target_col = "Werkelijkheid",
    prediction_col = "Prediction",
    counts_col = "n",
    palette = "Blues",
    add_sums = TRUE,
    theme_fn = ggplot2::theme_light,
    sums_settings = sum_tile_settings(
      palette = "Greens",
      label = "Totaal",
      tc_tile_border_color = "black"
    )) +
    
    # Customize the labels
    labs(
      x = "Werkelijke uitkost",
      y = "Voorspelde uitkomst",
      caption = sCaption
    ) +
    
    Set_LTA_Theme()
  
  # Add LTA elements
  confusion_plot <- Add_LTA_Theme_Elements(confusion_plot, 
                                           title_subtitle = TRUE)
  
  print(confusion_plot)
Figuur 4.11: Confusion matrix ten opzichte van Retentie na 1 jaar

4.6.4 Uitleggen of verklaren?

Naast de accuraatheid van het model is het ook belangrijk om te weten welke factoren het meest bijdragen aan de voorspelling van retentie na 1 jaar. Daarin gaat de vergelijking met de prestaties van het basemodel mank. Dat model geeft op geen enkele manier aan waarom een student een kans op succes heeft, anders dan - ‘dit is gebruikelijk in deze opleiding’.

Ongeacht de mate van accuraatheid, is het voor onderzoek naar kansengelijkheid essentieel om te weten welke factoren het meest bijdragen aan de voorspelling van retentie na 1 jaar. Het gaat erom dat we het belang van de factoren in de voorspellingen kunnen begrijpen en duiden. Machine Learning is hiervoor uitstekend geschikt, omdat het de mogelijkheid biedt om de belangrijkste factoren en hun invloed te leren kennen (Shmueli, 2010; Shmueli & Koppius, 2011).

 

Verantwoording

Deze analyse maakt deel uit van het onderzoek naar kansengelijkheid van het lectoraat Learning Technology & Analytics van De Haagse Hogeschool: No Fairness without Awareness | Het rapport is door het lectoraat ontwikkeld in Quarto 1.6.39. | Template versie:

 

Copyright

Dr. Theo Bakker, Lectoraat Learning Technology & Analytics, De Haagse Hogeschool © 2023-2025. Alle rechten voorbehouden.